#main script using package findmfdia
rm(list=ls())
detach("package:findMFDIA")
library(findMFDIA)

source("../R/findpeakgroups.R")
source("../R/plotting.R")
librarydir<-"~/prog/findmfdia/data"
datadir <- "/media/data/testdata/AquaWater"
datadir <- "/media/data/testdata/AquaHuman"
datadir <- "/media/data/testdata/AQUAYest"



dir(datadir)
files <- dir(datadir,"*sqlite")
length(files)

### load irts
irttrans2 = getTransitions(file.path(librarydir,"DIA_iRT.tsv"))
irttrans2 = filterProdByMass(irttrans2 , 200)
## load test
alltrans=getTransitions(file.path(librarydir,"AQUASky_ShotgunLibrary_3t_345_sh.10corrected.csv"))
alltrans = filterProdByMass(alltrans , 200)

decoys <- NULL
for(fcount in 1:30)
    {
        fcount <- 10
        fname <- files[fcount]
        con=createConnection(file.path(datadir,fname))
        
        tmp <- strsplit(fname,"_")[[1]]
        dataset <- tmp[2]
        dataset
        dil <- as.numeric(tmp[3])
        dil
        
        addIndexMZRT2Features(con)
        ## finding peak groups ######
        createBestFeaturesView(con,4,2,200,"ms1view")
        createBestFeaturesView(con,0,0,5)

        newtrans <- calibrateRTTransitions(con,irttrans2,alltrans,rterror=3,relax=1)
        
        ## plots new against old
        x <- unlist(lapply(alltrans,function(x){x$rt}))
        y <- unlist(lapply(newtrans,function(x){x$rt}))
        plot(x,y)
        
        createBestFeaturesView(con,4,2,50)
        scorex <- findTransitions(newtrans,con,mzerror=0.05,extRTerr=180,rterror=8,relax=0)

        length(scorex$group)

        if(length(scorex)==0)
            next
        par(mfrow=c(1,1))
        score <- scorex$score

        length(scorex$group[score$decoy==0])
        length(scorex$group[score$decoy==1])
        
        visTransitions(con,scorex$group[score$decoy==0],sleep=0.2,main="target")
        #visTransitions(con,scorex$group[score$decoy==1],sleep=1.2,main="decoy")
        
        tsc <- score[score$decoy==0,]
        groups <- scorex$groups[score$decoy==0]
        dsc <- score[score$decoy==1,,drop=FALSE]
        cat("dil ", dil ,"targets ",dim(tsc)[1],"decoys ",dim(dsc)[1],"\n")
        decoys <- rbind(decoys,c(fname,dim(dsc)[1]))

        ## visulize duplicates
        sc <- tsc
        d<-sc$sequence[duplicated(sc$sequence)]
        

        if(1){
            for(i in d){
                x <- which(sc$sequence==i)
                print(scorex$score[x,])
                par(mfrow=c(1,2))
                plotGroup(con,groups[[x[1]]]$trans,groups[[x[1]]]$featuregroups)
                plotGroup(con,groups[[x[2]]]$trans,groups[[x[2]]]$featuregroups)
                #Sys.sleep(3)
            }
        }
        x <- which(duplicated(sc$sequence))
        ##remove duplicates from dataset
        if(length(x)>0){
            utsc <- tsc[-x,]
        }else{
            utsc <- tsc
        }
        dim(utsc)
        cat("dil",dil,"unique targets ",dim(utsc)[1],"\n")
       
        gg <- toSimpleFormat(utsc)
        stemp <- strsplit(fname,"\\.")[[1]][1]
        write.table(gg,file=file.path(datadir,paste(stemp,".txt",sep="")))

        res <- tryCatch(
            dbClearResult(dbListResults(con)[[1]])
            )
        res <- tryCatch(dbDisconnect(con))
    }


#x <- "/path/to/file.txt"
#y<-strsplit(x,"/")
#y
#y <- y[[1]]
#y <- y[length(y)]
#strsplit(y,"\\.")[[1]]

#file.stem <- function(x){
#        bn <- basename(x)
#        gsub("\\..*$", "", bn)
#}
#file.stem("/the/path/to/afile.txt")
